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ABSTRACT 

This  report  presents  a  set  of  test  problems  for  nonlinear  equations  and  non¬ 
linear  least-squares  algorithms.  These  problems,  sent  to  us  by  C.V.  Nelson  of  the 
Maine  Medical  Center,  come  from  a  dipole  model  of  the  heart.  They  are  6x6  or 
8x8,  easy  to  code,  cheap  to  evaluate,  and  not  easy  to  solve.  In  support  of  the 
latter  contention,  we  present  test  results  from  MINPACK  and  NL2SOL. 


1.  Introduction 

Consider  the  following  two  problems: 

NLEQ:  Given  F  :Rn-RB,  solve  F(x)  =  0. 

NLLS:  Given  F:RB-Rm,  minimize  <K*)  =  y  F(x)TF(x). 

Most  algorithms  for  the  solution  of  these  two  problems  are  based  on  the  assumption  that  F 
can  be  adequately  modeled  by  an  affine  function  in  some  neighborhood  of  a  point  of  interest^ 
whether  that  point  is  close  to  or  far  away  from  the  solution  to  the  problem.  The  purpose  of  this 
note  is  first  to  add  an  interesting  test  function  to  the  current  list  of  test  problems  for  nonlinear 
equations  and  nonlinear  least  squares,  and  second  to  use  that  test  function  to  give  some  indications 
that  affine  modeling  by  the  first  two  terms  of  the  Taylor  series  is  not  necessarily  the  best  strategy 
for  Newton-type  methods  far  from  the  solution  to  the  problem.  This  claim  will  be  supported  by 
computational  results  from  five  codes: 

•  HYBRD  from  MINPACK  —  intended  for  NLEQ  —  uses  a  secant  affine  model  for  F; 

•  f-d  HYBRD  from  MINPACK  —  intended  for  NLEQ  —  uses  a  Newton  affine  model  for  F; 

•  LMDIF  from  MINPACK  —  intended  for  NLLS  —  uses  a  Newton  affine  model  for  F; 

•  NL2SNO  from  NL2SOL  —  intended  for  NLLS  —  uses  a  mixed  quadratic  model  for  <Jr, 

•  DN2F  —  intended  for  NLLS  —  a  slight  modification  of  NL2SOL. 

In  the  next  section,  we  will  give  a  brief  description  of  the  test  problems  mentioned  above. 
Section  3  will  discuss  the  way  the  five  codes  handle  their  corresponding  problems  and  Section  4 
gives  the  computational  results  together  with  some  conclusions.  A  listing  of  the  Fortran  subroutine 
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used  to  evaluate  the  test  function  appears  in  the  appendix,  as  do  five  sets  of  experimental  data  and 
corresponding  solutions. 

2.  The  Test  Problem 

We  give  here  a  brief  description  of  a  problem  communicated  to  us  by  C.V.  Nelson  (1980). 
The  interested  reader  can  find  a  detailed  description  of  the  problem  together  with  the  derivation  of 
the  equations  in  Nelson  and  Hodgkin  (1981).  The  problem  comes  from  experiments  using  two 
artificial  dipoles  in  a  circular  disk  containing  electrolyte  to  determine  the  resultant  dipole  moment 
of  the  human  heart.  Potentials  from  electrodes  around  the  boundary  of  the  disk  are  measured  and 
the  Gabor-Nelson  (1954)  equations  are  used  to  solve  for  the  magnitude,  directions,  and  locations 
of  the  two  independent  dipoles  in  the  disks.  The  problem  reduces  to  that  of  solving  the  following 
8x8  system  of  nonlinear  equations  in  the  unknowns  a,  b,  c,  d,  t,  u,  v,  w : 

a  +  b  =  2  Mx 

c  +  d  —  2  My 

ta  +  ub  -  vc  -  wd  =  2  A 

va  +  wb  +  tc  +  ud  =  2  B 

ait1—  v2)  —  2crv  +  b(u2—  w2)  -  2 duw  =  2  C 

c(r2-  v2)  +  2 atv  +  d{u2-  w2)  +  2buw  =  2  D 

at(t 2-  3v2)  +  cv(v2-  3f2)  +  bu(u 2-  3vv2)  +  dw{w2-  3 u2)  =  IE 

ct{t2-  3v2)  -  av(v2-  3r2)  +  du{u2-  3w2)  -  bw(w2-  3 u2)  =  IF, 

where  the  right-hand  sides  of  the  equations  in  each  experiment  were  evaluated  from  measured 
potentials.  Note  that  since  the  first  two  equations  are  linear  in  the  variables,  we  can  eliminate  two 
of  the  first  four  variables  to  obtain  a  6x6  system  which  might  be  easier  to  solve.  Below,  we  will' 
assume  that  the  equations  have  been  rewritten  to  have  right  hand  side  equal  to  zero  and  we  will 
refer  to  the  8X8  system  as  the  full  problem  and  to  the  6x6  system  as  the  reduced  problem. 

This  appears  to  be  a  good  test  problem  for  NLEQ  and  NLLS  algorithms  because  it  is  easy  to 
code,  cheap  to  evaluate,  and  hard  to  solve,  judging  from  our  experience  with  the  five  codes  listed 
above.  In  all  of  the  problems,  we  did  a  LINPACK-type  condition  estimate  of  the  Jacobian  matrix 
F'(x)  at  the  solution  and  it  was  never  greater  that  105. 

3.  The  Codes  Tested 

In  this  section  we  will  give  a  brief  description  of  the  five  codes  we  tested.  We  used  the 
finite-difference  version  of  each.  All  five  codes  use  model-trust  region  algorithms;  for  detailed 
explanations  of  such  algorithms,  see  Chapters  6,  8,  and  10  of  Dennis  and  Schnabel  (1983). 

Only  two  of  the  five  codes  tested,  HYBRD  and  f-d  HYBRD,  use  the  nXn  structure  of  the 
problem.  They  both  use  a  modification  of  Powell’s  (1970)  dogleg  algorithm  as  a  global  strategy, 
and  the  local  strategy  is  based  on  an  affine  model  of  F{x).  In  f-d  HYBRD  this  is  done  using  for¬ 
ward  differences.  In  HYBRD,  the  model  is  constructed  at  each  iterate  using  the  Broyden  (1965) 
secant  update  to  approximate  the  Jacobian  F'(x).  Broyden’s  scheme  provides  a  model  that  matches 
the  two  most  current  F-values  rather  than  the  current  F  and  F '  values.  This  indicates  that  perhaps 
the  Broyden  scheme  is  better  able  to  remember  the  general  shape  of  the  function  over  several  itera¬ 
tions  far  from  the  solution,  and  conversely,  less  able  to  forget  outdated  information  as  the  iteration 
proceeds.  The  initial  model  Jacobian  is  provided  by  finite  differences  and  if  the  updated  affine 
model  seems  to  be  too  inaccurate  for  reasonable  progress  at.  any  iteration,  then  the  Jacobian 
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approximation  is  refreshed  by  finite  differences.  When  the  Broyden  update  is  used,  the  linear  alge¬ 
bra  to  compute  an  iterate  only  requires  0(n2)  flops. 

The  other  three  programs  consider  the  nXn  systems  as  nonlinear  least-squares  problems,  i.e., 
instead  of  solving  F(x )  =  0,  they  solve  min  y F(x)TF(x ).  By  posing  the  problems  in  this  way,  we 

are  able  to  use  the  LMDIF  code  from  MINPACK  and  the  NL2SNO  code  from  the  NL2SOL  pack¬ 
age  to  obtain  a  solution.  As  a  result  of  these  experiments,  we  discovered  that  NL2SOL  was  com¬ 
puting  a  lower  bound  for  the  Marquardt  parameter  that  was  smaller  than  intended.  DN2F  is  a  vari¬ 
ant  of  NL2SOL  that  appears  in  PORT  3,  the  third  edition  of  the  subroutine  library  described  by 
Fox,  Hall,  and  Schryer  (1978).  It  is  similar  to  NL2SOL,  except  that  the  bounds  on  the  Marquardt 
parameter  are  computed  as  in  LMDIF,  the  trust  radius  is  sometimes  not  changed  at  the  beginning 
of  an  iteration  (specifically,  the  trust  radius  is  not  allowed  to  decrease  if  the  last  step  was  a  “good” 
one),  and  scaling  (i.e.,  the  choice  of  Dk  in  (3.1))  is  based  on  the  infinity-norm  of  the  Jacobian 
matrix  columns  rather  than  their  Euclidean  norm.  AD  three  of  these  codes  require  0(n 3)  flops  to 
do  the  linear  algebra  at  each  iteration.  LMDIF  uses  a  modification  of  the  Levenberg-Marquardt 
algorithm,  i.e.,  given  a  current  estimate  xk  to  the  solution,  it  determines  a  search  direction  s k  and 
subsequent  iterate  xk+l  by 

U\Jk  +  MDk)sk  =  -J\F(xk)  (3.1) 

**  +  l  =  xk  +  sk 

where  Jk=F'(xk),  Dk  is  a  positive  diagonal  matrix,  and  p,*  is  a  nonnegative  scalar  that  is  adap¬ 
tively  chosen. 

Thus  LMDIF  builds  a  local  quadratic  model  of  (j>  by  matching  the  current  functional  and  gra¬ 
dient  values  and  using  J(xk)TJ(xk )  as  a  Hessian  approximation.  We  call  this  the  Gauss-Newton 
model  of  4>,  and  it  is  the  quadratic  model  obtained  by  building  a  Newton  affine  model  of  F(x)  at 
xk,  i.e.,  a  model  that  matches  F  and  F'  at  xk,  and  then  taking  the  sum  of  squares  of  the  affine 
model.  This  Gauss-Newton  model  of  d>  only  coincides  with  the  fuU  Newton  quadratic  model  of  <$> 
if  F  is  truly  affine  or  zero  in  each  residual.  However,  the  problems  being  considered  here  have 
zero  residuals  at  the  solution,  so  the  Gauss-Newton  method  wiU  have  the  same  quadratic  conver¬ 
gence  as  Newton's  method  near  the  solution. 

The  Hessian  approximation  used  by  NL2SNO  and  DN2F  includes  a  cheap  variable-metrit' 
secant  approximation  to  the  part  of  the  Hessian  of  (J>  that  the  Gauss-Newton  model  neglects.  Thus, 
it  is  a  compromise  between  the  Gauss-Newton  use  of  information  only  at  xk  and  the  secant 
method’s  memory  of  past  F-values.  The  algorithm  decides  adaptively  at  each  iteration  whether  to 
use  this  Hessian  augmentation.  The  decision  is  based  on  quadratic  model  information  from  a  tfust- 
region  implementation.  For  details  the  reader  should  consult  Dennis,  Gay,  and  Welsch  (1981)  or 
Dennis  and  Schnabel  (1983). 

Both  the  NLLS  algorithms  use  basically  the  same  global  strategy.  It  is  usually  viewed  as 
obtaining  the  next  iterate  by  minimizing  a  local  quadratic  model  of  <J)  in  some  region  about  the 
current  iterate  where  the  model  can  be  reasonably  trusted  to  adequately  model  <j>(*)-  The  only 
difference  is  that  LMDIF  never  has  to  worry  about  encountering  negative  curvature  in  the  model 
and  so  the  computation  of  each  iterate  is  a  bit  simpler.  Details  for  LMDIF  can  be  found  in  Mord 
(1978)  or  in  Chapter  6  of  Dennis  and  Schnabel  (1983).  Details  for  NL2SNO  can  be  found  in  Gay 
(1981). 

AU  three  implementations  compute  an  approximate  Jacobian  matrix  by  forward  differences 
using  a  step  size  in  the  /'th  coordinate  direction  of  \xi\(machep)x'2 ,  where  machep  is  the  machine 
epsilon  or  unit  rounding  error  of  the  arithmetic.  Also,  all  three  codes  have  the  capability  of  either 
letting  the  user  choose  a  fixed  diagonal  scaling-matrix  for  the  independent  variable,  or  determining 
it  internally  and  updating  it  at  each  iteration. 
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4.  Computational  Results 

Everything  is  coded  in  Fortran  and  run  on  a  VAX®ll/750  under  a  UNIX® operating  system 
using  DOUBLE  PRECISION  arithmetic  (for  which  we  take  machep  =  2-55).  Both  the  full  and 
reduced  problems  are  run  with  and  without  scaling  of  the  independent  variable  (i.e.,  with  Dk  in 
(3.1)  chosen  dynamically  by  the  algorithm  and  with  Dk  =  I).  Also,  we  execute  all  five  algorithms 
using  as  initial  guess  x0,  10x0.  and  10(k0,  where  x0  is  the  experimental  data  provided  by  Nelson 
and  given  in  the  appendix.  Since  there  are  five  sets  of  data,  this  gives  us  a  total  of  60  test  cases. 
We  always  found  the  same  single  solution  for  each  experiment,  independent  of  initial  guess,  scal¬ 
ing  option,  and  whether  or  not  we  solved  the  full  or  reduced  system.  In  comparing  off-the-shelf 
codes,  there  is  always  the  problem  of  differing  stopping  conditions.  We  used  the  same  tolerances 
for  the  same  tests  whenever  possible.  The  convergence  tests  used  by  the  NLLS  codes  are  quite 
similar,  but  our  version  of  HYBRD  uses  only  a  test  on  the  relative  change  in  the  dependent  vari¬ 
able.  The  other  codes  use  this  test  also;  the  tolerance  was  set  to  machep 1/2  for  all  the  runs. 

To  obtain  the  results  reported  for  HYBRD,  we  had  to  comment  out  its  tests  for  lack  of  pro¬ 
gress.  Without  this  change,  HYBRD  failed  on  about  a  third  of  the  test  problems.  We  did  not 
have  to  make  this  change,  however,  to  f-d  HYBRD,  which  we  had  obtained  by  modifying  HYBRD 
to  force  a  new  finite-difference  calculation  of  the  Jacobian  after  every  successful  step. 

Output  from  NL2SNO  and  DN2F  tells  us  that  the  algorithms  tended  to  use  the  augmented 
Gauss-Newton  model  far  away  from  the  solution  and  the  traditional  Gauss-Newton  model  near  the 
solution.  This  seems  to  suggest  that  the  augmentation  of  the  Gauss-Newton  Hessian  not  only  plays 
an  important  role  in  large-residual  problems,  but  also  gives  a  better  approximation  to  the  least- 
squares  Hessian  in  the  small-residual  case  when  we  are  far  from  the  solution.  Perhaps  our  earlier 
discussion  of  the  possible  advantages  of  memory  in  the  secant  methods  is  relevant  here. 

On  the  other  hand,  the  occasional  efficiency  of  DN2F  suggests  that  sometimes  it  would  be 
worthwhile  in  the  nonlinear  equations  problem  to  find  an  analog  to  the  NL2SOL  secant  augmenta¬ 
tion,  sizing,  and  model  switching  so  that  memory  of  past  points  can  be  suppressed  as  the  solution  is 
approached,  while  providing  a  better  model  than  the  affine  approximation  for  the  function  in  the 
early  iterations.  A  sophisticated  new  approach  to  this  problem  based  on  tensor  updating  is  sug¬ 
gested  in  Schnabel  and  Frank  (1984). 
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Table  1  -  Total  number  of  residual  calculations  with  scaling. 


HYBRD 

LMDIF 

NL2SNO(*) 

DN2F(*) 

f-d  HYBRD 

Full 

Prob. 

Red. 

Prob. 

Full 

Prob. 

Red. 

Prob. 

Full 

Prob. 

Red. 

Prob. 

g] 

Full 

Prob. 

Red. 

Prob. 

*0 

123 

327 

217 

286 

H 

I 

10*0 

188 

211 

381 

170 

346 

186 

337 

192 

100*0 

170 

104 

361 

219 

391 

228 

223 

325 

196 

*0 

37 

29 

45 

35 

37 

29 

II 

10*o 

1292 

61 

108 

91 

109 

m 

100*o 

1902 

356 

4302 

1407 

111 

*0 

311 

273 

1348 

696 

1454 

705 

295 

448 

m 

10*o 

2452 

863 

3525 

963 

4904 

1108 

4126 

610 

428 

419 

lOQxo 

4278 

1916 

5338 

3684 

5035 

7117 

451 

577 

■ 

n 

626 

1659 

569 

1303 

551 

211 

257 

D 

10*o 

795 

661 

755 

3063 

2723 

IB 

E9 

i 

E9 

1183 

5199 

809 

302 

v  305  " 

*0 

385 

419 

508 

1747 

n 

545 

1289 

V 

10*0 

IB 

978 

630 

1082 

459 

178 

228 

100*o 

232' 

a 

H 

226 

291 

I:  Exp.  791129  III:  Exp.  0121a  V:  Exp.  0121c 

H:  Exp.  793226  IV:  Exp. 0121b 

(*)  Both  NL2SN0  and  DN2F  were  run  with  the  LMDIF  default  initial  step  bound  of 
100*  ||Z>o*Jtoll-  Also  V(RDFCMX),  the  maximum  factor  by  which  the  trust  region  radius 
may  be  increased  at  one  time  is  changed  from  the  default  value  of  4  to  be  the  same  as  the 
LMDIF  default  value  of  2. 


January  31,  1986 


6- 


Table  2  -  Total  number  of  residual  calculations  without  scaling. 


HYBRD 

LMDIF 

NL2SNO(*) 

DN2F(*) 

f-d  HYBRD 

■ 

1 

Red. 

Prob. 

Red. 

Prob. 

H 

M 

1 

1 

1 

Red. 

Prob. 

*0 

180 

B 

B 

443 

552 

344 

j^jjj 

mn 

262 

I 

10*0 

195 

EB 

— 

205 

1269 

183 

274 

172 

289 

207 

100*Q 

144 

a 

219 

D 

309 

222 

327 

196 

■ 

■ 

37 

29 

45 

29  ’ 

109 

138 

1115 

356 

D 

85 

■ 

401 

2166 

683 

9027 

1594 

664 

335 

523 

145 

*0 

1904 

500 

B 

m 

3924 

2039 

4144 

2288 

397 

727 

100*  0 

716 

13736 

1990 

4818 

7116 

H 

399 

B 

B 

1228 

506 

1284 

419 

461 

El 

El 

E9 

m 

204 

403 

100*o 

870 

B 

B 

815 

B 

■ 

B 

B 

B 

529 

2153 

905 

1264 

445 

163 

151 

9 

10*o 

1006 

375 

939 

1041 

4596 

m 

840 

834 

158 

2083 

1 

B 

wm 

B 

23801 

2749 

2320 

B 

144 

I:  Exp.  791129  ffl:  Exp.  0121a  V:  Exp.  0121c 

II:  Exp.  791226  IV:  Exp.  0121b 

(*)  Both  NL2SN0  and  DN2F  were  run  with  the  LMDIF  default  initial  step  bound  of  100*  ||*0||. 

Also  V(RDFCMX),  the  maximum  factor  by  which  the  trust  region  radius  may  be  increased 
at  one  time  is  changed  from  the  default  value  of  4  to  be  the  same  as  the  LMDIF  default 
value  of  2. 
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Let2  =  (2M„  ^My,  2  A,  ££,  £  C,  2  D ,  2  E,  2 F)T • 

Exp.  791129 

2  =  (.485,  -.0019,  -.0581,  .015,  .105,  .0406,  .167,  -.399)r 
X0  =  (.299,  .186,  -.0273,  .0254,  -  .474,  .474,  -.0892,  .0892)7 

X *  ~  (  —  6.321349025c  —  3,  4.913213490c-l,  -  1.998156408e -3,  9. 8 15640840c -5 

1.226569755c -1,  - 1. 003153205c  -  1,  -4.023517593e +0,  -2.071785527e -2)r 


Exp.  791226 

2  =  (-.69,  -.044,  -1.57,  -1.31,  -2.65,  2.0,  -12.6,  9.48)r 
X0  =  (-.3,  -.39,  .3,  -  .344,  -1.2,  2.69,  1.59,  -1.5)r 

X *  =  (  —  3.116266056c  —  1,  -3.783733944c  - 1,  3.282442301c -1,  -3.722442301e  - 1 
-  1.282227094c +0,  2.494300312e  +  0,  1.554865879c +0,  -  1.384637843c +0)r 


Exp.  0121a 

2  =  (-.816,  -.017,  -1.826,  -.754,  -4.839,  -3.259,  -14.023,  15.467)T 
X0  =  (-.041,  -.775,  .03,  -.047,  -2.565,  2.565,  -.754,  .754)T 

X*  ~  (3.099869097c -3,  -8.190998691c-l,  -2.239405352c -4,  -  1.677605946c -2 

2.681514498e  +0,  2.250215931c +0,  -2.024170463c +  1,  7.970982952c -l)r  ; 


Exp.  0121b 

2  =  (-.809,  -.021,  -2.04,  -.614,  -6.903,  -2.934,  -26.328,  18.639)r 
X0  =  (-.056,  -.753,  .026,  -.047,  -2.991,  2.991,  -.568,  .568)r 

X*  ~  (9.034542990e  -3,  -8.180345430e  - 1,  -4.450738446c -4,  -2.055492616c -2 
2. 773429036c  +0,  2. 529477259c +0,  - 1 .480097186c  + 1,  5. 220468844c  -1)T 

Exp.  0121c 

2  =  (-.807,  -.021,  -2.379,  -.364,  -10.541,  -1.961,  -51.551,  21.053)7 
X0  =  (-.074,  -.733,  .013,  -  .034,  -3.632,  3.632,  -  .289,  .289)r 

X*  ~  (5.140417418e -2,  -8.584041742c - 1 ,  1. 047333626c -3,  -2.204733363c -2 
2.861205288e  +0,  2.949155438e +0,  - 8. 304243489c +0,  - 1.4549924 13e-l)T 
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SUBROUTINE  FCN(N,  X,  FVEC,  IFLAG) 

C 

C  This  is  an  example  of  a  subroutine  to  evaluate  the  "heart  model" 
function  to  use  with  HYBRD  of  MINPACK. 

IPICK  ■  1  - >  full  problem. 

IPICK  *  2  - >  reduced  problem. 

INTEGER  N,  IFLAG 
DOUBLE  PRECISION  X(N) ,FVEC(N) 

C 

INTEGER  IPICK 

DOUBLE  PRECISION  SIGMAX,  SIGMAY,  SIGMAA,  SIGMAB,  SIGMAC,  SIGMAD, 

1  SIGMAE,  SIGMAF 

COMMON  /SIGMA/  SIGMAX,  SIGMAY,  SIGMAA,  SIGMAB,  SIGMAC,  SIGMAD, 

1  SIGMAE,  SIGMAF,  IPICK 

C 

DOUBLE  PRECISION  AT,  ATV ,  AV,  B,  BU,  BUW,  BW,  CT,  CTV,  CV,  D,  DU, 

1  DUW,  DW,  T2M3V2,  T2MV2,  U2M3W2,  U2MW2 ,  V2M3T2,  W2M3U2 

C 

IF  (IPICK  .NE.  1)  GO  TO  10 
C 

T2MV2  «  X(5 ) **2  -  X(7)**2 

U2MW2  *  X(6)**2  -  X( 8 ) **2 

T2M3V2  «  X(5 ) **2  -  3 . 0D0  *  X(7)**2 

V2M3T2  *  X( 7 ) **2  -  3 . 0D0  *  X(5)**2 

U2M3W2  «  X(6) **2  -  3 . 0D0  *  X<8)**2 

W2M3U2  ■  X ( 8 ) **2  -  3.0D0  *  X(6)**2 

CTV  *  X { 3 V  *  X{ 5 )  *  X(7) 

DUW  *  X( 4 )  *  X ( 6 )  *  X ( 8 ) 

ATV  «  X( 1 )  *  X ( 5 )  *  X(  7 ) 

BUW  b  X( 2 )  *  X(6)  *  X{ 8) 

AT  *  X( 1)  *  X(5) 

BU  =  X( 2 )  *  X( 6 ) 

CV  *  X(3)  *  X ( 7 ) 

DW  *  X(  4  )  *  X( 8 ) 

CT  *  X(  3  )  *  X(  5  ) 

AV  «  X( 1)  *  X{7) 

DU  *  X ( 4 )  *  X { 6 ) 

BW  ■  X( 2 )  *  X ( 8 ) 

C 

FVEC ( 1 )  *  X(1)  +  X ( 2 )  -  SIGMAX 
FVEC ( 2 )  »  X( 3 )  +  X(4 )  -  SIGMAY 
FVEC ( 3 )  b  AT  +  BU  -  CV  -  DW  -  SIGMAA 
FVEC ( 4 )  b  AV  +  BW  +  CT  +  DU  -  SIGMAB 

FVEC ( 5 )  b  X{ 1 ) *T2MV2  -  2.0DO*CTV  +  X(2)*U2MW2  -  2.0D0*DUW  -  SIGMAC 

FVEC ( 6 )  =  X( 3 ) *T2MV2  +  2.0D0*ATV  +  X(4)*U2MW2  +  2.0D0*BUW  -  SIGMAD 

FVEC ( 7 )  »  AT*T2M3V2  +  CV*V2M3T2  +  BU*U2M3W2  +  DW*W2M3U2  -  SIGMAE 

FVEC { 8 )  b  CT*T2M3V2  -  AV*V2M3T2  +  DU*U2M3W2  -  BW*W2M3U2  -  SIGMAF 

GO  TO  999 
C 

TO  T2MV2  ■  X { 3 ) **2  -  X{5)**2 
U2MW2  s  X(4 ) »*2  -  X ( 6 ) **2 
T2M3V2  b  X ( 3 ) **2  -  3 . ODO  *  X(5)**2 
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V2M3T2  ■  X ( 5 ) * *2  -  3.0D0  *  X<3)**2 
U2M3W2  ■  X(4)«*2  -  3 . 0D0  *  X(6)**2 
W2M3U2  ■  X{6)*#2  -  3.0D0  *  X{4)**2 
B  ■  SIGMAX  -  X(1) 

D  •  SIGMAY  -  X(2) 


CTV 

m 

X ( 2 )  *  X{3)  *  X( 5) 

DUW 

m 

D  *  X{ 4  )  *  X{6 ) 

ATV 

m 

X( 1 )  *  X(3)  *  X(5) 

BUW 

m 

B  *  X ( 4 )  *  X{  6  ) 

AT 

u 

X(1)  *  X ( 3 ) 

BU 

m 

B  •  X( 4 ) 

CV 

* 

X ( 2 )  *  X(5) 

DW 

m 

D  *  X(  6 ) 

CT 

m 

X( 2 )  *  X( 3 ) 

AV 

* 

X(1)  *  X ( 5 ) 

DU 

ft 

D  *  X  {  4  ) 

BW 

ft 

B  X(  6  ) 

C 

FVEC ( 1 )  *  AT  +  BU  -  CV  -  DW  -  SIGMAA 

FVEC ( 2 )  ■  AV  +  BW  +  CT  +  DU  -  SIGMAB 

FVEC { 3 )  «  X(1)*T2MV2  -  2.0DO*CTV  +  B*U2MW2  -  2.0D0*DUW  -  SIGMAC 

FVEC ( 4 )  -  X{ 2 ) *T2MV2  +  2.0D0*ATV  +  D*U2MW2  +  2.0D0*BUW  -  SIGMAD 

FVEC ( 5 )  *  AT*T2M3V2  +  CV*V2M3T2  +  BU*U2M3W2  +  DW*W2M3U2  -  SIGMAE 
FVEC ( 6 )  «  CT*T2M3V2  -  AV*V2M3T2  +  DU*U2M3W2  -  BW*W2M3U2  -  SIGMAF 
999  RETURN 
C 

END 
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